clear
rng default
addpath /Users/tew647/mosek/9.3/toolbox/r2015aom
%% Load useful matrices
load PXY_cond 
PYX_cond_split=[PXY_cond(1:3); PXY_cond(4:end)];
load V_CS

%% Construct grid of parameters to check for each realisation of X
gran=0.01;
u0grid_temp=0; %[Xchosen=0, Y=1] [Xchosen=0, Y=2], location normalisation 
u1grid_temp_all=[V_11_CS V_12_CS]; %[Xchosen=1, Y=1] [Xchosen=1, Y=2], scale normalisation
u2grid_temp=-20:gran:20;  %[Xchosen=2, Y=1] [Xchosen=2, Y=2]
% try different values of gran and different sizes of grid hypercube as discussed in Appendix C of the paper (need powerful cluster to run the code for higher values of gran and larger sizes of hypercube)

%% Construct identified set for each realisation of Y
ntypes_w=2;
IdSetW=cell(ntypes_w,1);

for x=1:ntypes_w %we can proceed separately across types as long as we do not impose independence from Y
    u1grid_temp=u1grid_temp_all(x);
    [ca, cb, cc] = ndgrid(u0grid_temp, u1grid_temp, u2grid_temp);
    u0grid=ca(:);
    u1grid=cb(:);
    u2grid=cc(:);
    sg=size(u0grid,1);
    % Components invariant across parameter values
    [sgb, ...
    equalorder, ...
    Ugridreduced, sgreduced,...
    id_Ugridreduced,...
    s_id_gridreduced, sU_id_gridreduced,...
    n_U_sets, Tcoord, indices_pairs,...
    numeqconstraints,numzeros,coordrowseq,d3, e5, f5, g3, h5, i5, l3,m5, n5, fillAeq]=useful_anyparam3(u0grid, u1grid, u2grid, sg);  
    % Loop over parameter values, order of columns: [u0x u1x u2x]
    IdSetW{x}=OneSide3(PYX_cond_split(x,:), u0grid, u1grid, u2grid,...
                        sgb, ...
                        equalorder, ...
                        Ugridreduced, sgreduced,...
                        id_Ugridreduced,...
                        s_id_gridreduced, sU_id_gridreduced,...
                        n_U_sets, Tcoord, indices_pairs,...
                        numeqconstraints,numzeros,coordrowseq,d3, e5, f5, g3, h5, i5, l3,m5, n5, fillAeq);  
end

[cx1, cx2] = ndgrid(1:size(IdSetW{1},1), 1:size(IdSetW{2},1));
cx1=cx1(:);
cx2=cx2(:);
IdSetW_final=[IdSetW{1}(cx1,:) IdSetW{2}(cx2,:)]; 
save('IdSetW_final.mat', 'IdSetW_final')


